home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / libs / vopl / glvopl.lha / glvopl / src / lsfit.c < prev    next >
Encoding:
C/C++ Source or Header  |  1989-07-27  |  4.7 KB  |  301 lines
  1. #include "vopl.h"
  2.  
  3. /*
  4.  * degree
  5.  *
  6.  *    set the order for a least-squares fit
  7.  */
  8. degree(ord)
  9.     int    ord;
  10. {
  11.     plotdev.degree = ord;
  12. }
  13.  
  14. /*
  15.  * avefit
  16.  * 
  17.  *    Draw a line y = average y[i] - used for degree 0
  18.  */
  19. void
  20. avefit(x, y, n)
  21.     float    x[], y[];
  22.     int    n;
  23. {
  24.     float    a;
  25.     int    i, na;
  26.     
  27.     a = 0.0;
  28.     na = 0;
  29.     for (i = plotdev.startind; i < n; i += plotdev.arrayind) {
  30.         na++;
  31.         a += WhatY(y[i]);
  32.     }
  33.  
  34.     if (na > 0)
  35.         a /= (float)na;
  36.  
  37.     move2(WhatX(plotdev.axes[XIND].min), a);
  38.     draw2(WhatY(plotdev.axes[XIND].max), a);
  39. }
  40.  
  41. /*
  42.  * llsfit
  43.  *
  44.  *    linear least square fit - used for degree 1
  45.  *
  46.  *    y = a + b.x
  47.  *
  48.  *    b =    n.SUM(xi.yi) - SUM(xi).SUM(yi)
  49.  *             -------------------------------
  50.  *           n.SUM(xi.xi) - SUM(xi).SUM(xi)
  51.  *
  52.  *    a =    yave - b.xave
  53.  */
  54. void
  55. llsfit(x, y, n)
  56.     float    x[], y[];
  57.     int    n;
  58. {
  59.     float    sumx, sumy, sumxy, sumxx, a, b, yave, xave, an;
  60.     int    i;
  61.     
  62.     sumx = sumy = sumxy = sumxx = 0.0;
  63.  
  64.     for (i = plotdev.startind; i < n; i += plotdev.arrayind) {
  65.         sumx += WhatX(x[i]);
  66.         sumy += WhatY(y[i]);
  67.         sumxy += WhatX(x[i]) * WhatY(y[i]);
  68.         sumxx += WhatX(x[i]) * WhatY(x[i]);
  69.     }
  70.  
  71.     b = (n * sumxy - sumx * sumy) / (n * sumxx - sumx * sumx);
  72.     yave = sumy / n;
  73.     xave = sumx / n;
  74.     a = yave - b * xave;
  75.  
  76.     yave = a + WhatX(plotdev.axes[XIND].min) * b;
  77.     move2(WhatX(plotdev.axes[XIND].min), yave);
  78.     yave = a + WhatX(plotdev.axes[XIND].max) * b;
  79.     draw2(WhatX(plotdev.axes[XIND].max), yave);
  80. }
  81.  
  82. static float     *b, *c, *d;
  83.  
  84. /*
  85.  *  orthofit
  86.  *
  87.  *   Calculate the orthogonal polynomial of the given degree 
  88.  *   through the data set (x, y) - used for degree 2 and above.
  89.  *
  90.  */
  91. int
  92. orthofit(xx, yy, degree, npnts)
  93.     float     *xx, *yy;
  94.     int     degree, npnts;
  95. {
  96.     
  97.     float    *x, *y, *pjm1, *pj;
  98.     float    ax, ay, xdiv; 
  99.     float     ortval();
  100.     int    j;
  101.     
  102.     /* get memory for variables */
  103.  
  104.     b = newm1(degree + 1);
  105.     c = newm1(degree + 1);
  106.     d = newm1(degree + 1);
  107.     pjm1 = newm1(npnts);
  108.     pj = newm1(npnts);
  109.     x = newm1(npnts);
  110.     y = newm1(npnts);
  111.  
  112.     for (j = 0; j < npnts; j++) {
  113.         x[j] = WhatX(xx[j]);
  114.         y[j] = WhatY(yy[j]);
  115.     }
  116.  
  117.     ortpol(npnts, x, y, pjm1, pj, degree);
  118.  
  119.     xdiv = (WhatX(plotdev.axes[XIND].max) - WhatX(plotdev.axes[XIND].min)) / (float)plotdev.precision;
  120.     ax = WhatX(plotdev.axes[XIND].min);
  121.  
  122.     ay = ortval(ax, degree);
  123.     move2(ax, ay);
  124.     for (j = 0; j <= plotdev.precision; j++) {
  125.          ax += xdiv;
  126.         ay = ortval(ax, degree);
  127.         draw2(ax, ay);
  128.     }
  129.     
  130.     free(b);
  131.     free(c);
  132.     free(d);
  133.     free(pjm1);
  134.     free(pj);
  135.     free(x);
  136.     free(y);
  137.  
  138.     return(j);
  139. }
  140.  
  141. /*
  142.  * ortpol
  143.  *
  144.  *    coeficients for orthogonal polynomial 
  145.  *
  146.  */
  147. int 
  148. ortpol(npnts, x, y, pjm1, pj, degree)
  149.     int    npnts, degree;
  150.     float    *x, *y, *pjm1, *pj;
  151. {
  152.     float    *s, *w;
  153.     float    p;
  154.     int    i, j;
  155.  
  156.     s = newm1(degree + 1);
  157.     w = newm1(npnts);
  158.  
  159.     for (j = 0; j <= degree; j++) {
  160.         b[j] = 0.0;
  161.         s[j] = 0.0;
  162.         d[j] = 0.0;
  163.     }
  164.  
  165.     for (i = 0; i < npnts; i++)
  166.         w[i] = 1.0;
  167.  
  168.     c[0] = 0.0;
  169.  
  170.     for (i = 0; i < npnts; i++) {
  171.         d[0] += y[i] * w[i];
  172.         b[0] += x[i] * w[i];
  173.         s[0] += w[i];
  174.     }
  175.  
  176.     d[0] /= s[0];
  177.  
  178.     if (degree == 1)
  179.         return(1);
  180.  
  181.     b[0] /= s[0];
  182.     for (i = 0; i < npnts; i++) {
  183.         pjm1[i] = 1.0;
  184.         pj[i] = x[i] - b[0];
  185.     }
  186.  
  187.     j = 0;
  188.  
  189.     for (;;) {
  190.         j++;
  191.         for(i = 0; i < npnts; i++) {
  192.             p = pj[i] * w[i];
  193.             d[j] +=  (y[i] - d[0]) * p;
  194.             p *= pj[i];
  195.             b[j] += x[i] * p;
  196.             s[j] += p;
  197.         }
  198.  
  199.             d[j] /= s[j];
  200.  
  201.         if (j == degree)
  202.             return(1);
  203.  
  204.         b[j] /= s[j];
  205.         c[j] = s[j] / s[j - 1];
  206.         for (i = 0; i < npnts; i++) {
  207.             p = pj[i];
  208.             pj[i] = (x[i] - b[j]) * pj[i] - c[j] * pjm1[i];
  209.             pjm1[i] = p;
  210.         }
  211.     }
  212. }
  213.  
  214. /*
  215.  * ortval
  216.  *
  217.  *    calculates the value of the othogonal polynomial
  218.  *    at the point x.
  219.  *
  220.  */
  221. float
  222. ortval(x, degree)
  223.     float     x;
  224.     int    degree;
  225. {
  226.     int     k;
  227.     float     prev, prev2,  val;
  228.     
  229.     k = degree;
  230.     val = d[k];
  231.     prev = 0.0;
  232.     while (k != 0) {
  233.         k--;
  234.         prev2 = prev;
  235.         prev = val;
  236.         val = d[k] + (x - b[k]) * prev - c[k + 1] * prev2;
  237.     }
  238.     return(val);
  239. }
  240.  
  241. /*
  242.  * calc_coef
  243.  *   
  244.  *     Note that the following two routines are not used to 
  245.  *     to interpolate for the polynomial fit of the function. They merely
  246.  *     to elicit the values of the coeficients of the polynomial.
  247.  *     
  248.  */
  249. calc_coef(degree, coef)
  250.     int    degree;
  251.     float    *coef;
  252. {
  253.         
  254.     float    **p;
  255.     int    i, j;
  256.  
  257.     p = newm2(degree + 1, degree + 1);
  258.  
  259.     for (i = 0; i <= degree; i++)
  260.         coef[i] = 0.0;
  261.  
  262.     p[0][0] = 1.0;
  263.     p[1][0] = -b[0];
  264.     p[1][1] = 1.0;
  265.  
  266.     coef[0] += d[0] * p[0][0] + d[1] * p[1][0];
  267.     coef[1] += d[1] * p[1][1];
  268.  
  269.     for (i = 2; i <= degree; i++) {
  270.         bp(i, b[i - 1], c[i - 1], p[i], p[i - 1], p[i - 2]);
  271.         for(j = 0; j <= i; j++) 
  272.             coef[j] += d[i] * p[i][j];
  273.     }
  274.  
  275.     fprintf(stderr, "coef:");
  276.     for (i = 0; i <= degree; i++)
  277.         fprintf(stderr, " %e\n", coef[i]);
  278.     fprintf(stderr, "\n");
  279.  
  280.     for (i = 0; i < degree + 1; i++)
  281.         free(p[i]);
  282.     
  283.     free(p);
  284. }
  285.  
  286. /*
  287.  * bp
  288.  *
  289.  */
  290. bp(degree, b, c, p3, p2, p1)
  291.     int     degree;
  292.     float    b, c;
  293.     float    *p1, *p2, *p3;
  294. {
  295.     int i;
  296.  
  297.     p3[0] = -b * p2[0] - c * p1[0];
  298.     for (i = 1; i <= degree; i++) 
  299.         p3[i] = p2[i - 1] - b * p2[i] - c * p1[i];
  300. }
  301.